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Abstract 

We show for a variety of classes of conservative PDEs that discrete gradient 
methods designed to have a conserved quantity (here called energy) also have 
a time-discrete conservation law. The discrete conservation law has the same 
conserved density as the continuous conservation law, while its flux is found by 
replacing all derivatives of the conserved density appearing in the continuous 
flux by discrete gradients. 

The discrete gradient method [3], [5l El [9] is a powerful and general technique for 
constructing integral-preserving integrators. It includes projection methods and the 
average vector field method [6] as special cases. In particular, it is useful for con- 
structing energy-preserving integrators. 

For PDEs, a conservation law is more fundamental than a conserved quantity, for 
the former implies the latter only in the presence of suitable boundary conditions, 
and the flux is (usually) a local quantity. Discrete gradient methods have been widely 
used to construct energy-preserving integrators for PDEs and their semidiscretiza- 
tions; see [2] for a survey. In this paper we show that discrete gradient methods also 
preserve energy conservation laws (ECLs) when the symplectic structure is constant. 
Discrete gradient methods are not usually symplectic or variational and the resulting 
conservation laws developed in this paper are not instances of Noether's theorem. 
Despite their nonvariational nature, the use of 'conservative' discretizations of sys- 
tems of conservation laws of mass, momentum and so on — that is, of discretizations 
with discrete conservation laws — is widespread in numerical methods for PDEs. The 
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typical case is that the conserved density is (one of) the dependent variable(s): for 
example, 

ut + f{u)^ = 

is a conservation law with conserved density u and flux f{u). Euler's method 

Ui - Uq 



At 



+ = 



is a time discretization in the form of a (time-) discrete conservation law. (The full 
discretization 

^ + A(/K)) 

where A is any difference operator in x, is a fully discrete conservation law, but the 
spatial discretization is covered by discrete versions of Noether's theorem (see, e.g., 
|1] and references therein) and is not the subject of this paper.) Another known 
case is when a conservation law arises in a system of Euler-Lagrange equations in 
which the Lagrangian admits a symmetry acting only on the dependent variables. 
If the symmetry is fairly simple, e.g., if it is linear, then there may exist symplectic 
or discrete Lagrangian methods that preserve the symmetry and which will have a 
discrete conservation law. See [8] where this situation is illustrated for the phase 
symmetry of the nonlinear Schrodinger equation. 

Discrete gradient methods are a way of ensuring that a time discretization has 
a conservation law when the conserved density, that we call energy in this paper, is 
not or cannot be one of the dependent variables, in a way that is independent of any 
variational structure or symmetry. 

A discrete gradient on a vector space V with inner product a^b is a function 
satisfying the discrete gradient axiom 

{z, - zoVVH{zo, z,) = H{z,) - H{zo) 

for all smooth functions H:V ^ and all Zq, Zi G V. For consistency in the 
continuous limit, generally one also requires 

VH{zo,zo) = VH{zo), 

although for conservation properties only the discrete gradient axiom is needed. A 
common discrete gradient is the Average Value discrete gradient given by 

VH{zo,z,)= [ Vif(e^i + (1 - O^o) de 
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In coordinates we write Hy for the ^/-component of VH. li z — {x,y), then the 
discrete gradient axiom becomes 

{xi - xofH^ + {yi - yofHy = H{xi, yi) - H{xo, yo). 

A discrete gradient method for the system Zt — K{z)VH{z), — —K, is a 
map Zq !->■ Zi given by 

Zi — Zq 



At 



K{zo,zi)VH(zo,zi) 



where = —K and (for consistency) K{zQ,Zf)) = K{zq,zq). 

We first give an exphcit calculation of the ECL for first-order canonical Hamilto- 
nian PDEs in one space dimension and their discrete gradient time discretizations. 

Proposition 1. Let z{x,t) e M", let K he an n x n antisymmetric matrix and let 
H be a Hamiltonian of the form H = J H{z, z^) <^x, and consider the equations of 
motion 

Zt^K—^K(H,-d,Hj. 
The discrete gradient method 



Zi — Zq 

At 



has a discrete ECL. Its conserved energy density, H , is the same as that of the PDF. 
Its flux is the same bilinear form as that of the continuous ECL, but with the gradients 
of the energy density replaced hy their discrete gradients. 

Proof. The continuous ECL can be found as follows. 

giving the energy conservation law + = where 

F = -HJKH^ + HjKd^H^ . 

Applying a discrete gradient method in {z, z^) to the PDE gives the time discretiza- 
tion 

^l^^KdT.-dJlZ). 



The change in energy density over one time step is 

H{zi, Zix) — H{zo, Zqx) -yt-'^ ^1 ~ , ■TJ^'^ ^Ix — Zqx 
: — tlz : r -H, : 

^ At At 



which estabhshes the result. □ 

Note that the discrete calculation exactly follows the continuous one, the key 
step being the discrete gradient axiom. The calculation is entirely local and does not 
require integration by parts. 



Example 1. The nonhnear wave equation = p, Pt = Qxx — V'{q) fits the frame- 
work of Proposition 1 with n = 2, K = and z = {q,p)^- The energy 
conservation law is 

H{q,p)t + F{q,p)x = Q 

where 

H = -p" + -ql + V{q), F = -pqx, 



and the discrete gradient ECL is 

H{qi,Pi) - H(qo,po) 
At 

where q ^ {q^ + qi)/2 and p = (po +Pi)/2. 



+ F{q,p)x = 



Proposition 2. Let K, be any constant (i.e., independent of z) skew-adjoint differ- 
ential operator and let % — H{z, Zx, Zxx, ■ ■ ■) be any differential function of z (and 
possibly X, although we suppress the x). Let K, and H determine the Poisson system 

zt = = ICE{H) 
oz 

where 

E{H) :— Hz — dxHz^ + dxxHz^^ — . . . 

is the Euler operator applied to H. The flux associated to the conserved density H 

is a quadratic form in the derivatives of H. A discrete gradient method applied to 
this system has a discrete ECL with conserved density H and flux given by the same 
quadratic form with derivatives of H replaced by discrete gradients. 
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Proof. From skew-adjointness of JC we can write 

JCb + JCa = dxS{a, h). 

(For example, S{a,b) = a^/C2& for /C = /Ci + K,2dx-) Then the flux for the energy- 
conservation law is the quadratic form [10] 



S{E{H), E{H)) - A{}CE{H), H) 

where 

oo n-l „ 
n=l 1=0 ^ 

The calculation of the discrete ECL now follows as in Proposition 1. □ 

In the literature there has been much attention paid to multisymplectic discretiza- 
tions of variational PDEs [T]. These often take the multi-Hamiltonian form 

Kzt + Lzx = VS'(z) 

which has a multisymplectic conservation law (d^; A Kdz)t + [dz A Ldz)^ = 0. Multi- 
symplectic methods preserve a discretization of this (differential) conservation law, 
but do not have a energy conservation law. (This is the PDE analogue of the fact 
that a integrator for a general Hamiltonian ODE cannot be both symplectic and 
energy- preserving) . Our next result covers this case. As in Proposition 2, it could be 
extended to differential operators K and arbitrary differential polynomials H, but 
the case considered here shows the basic structure. 

Proposition 3. Let % = H{z, z^) and let K he a (possibly singular) antisymmetric 
matrix. Then PDE 

dz 

has an ECL with conserved density H and flux given by a bilinear form in Zt and the 
derivatives of H . The discrete gradient method 

At z X 

has a discrete ECL with conserved density H and flux given by the same bilinear 
form with Zt replaced by {zi — ZQ)//S.t and the derivatives of H replaced by discrete 
gradients. 
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Proof. First we note that multi-Hamiltonian PDEs take the given form with H = 
S{z) — Lzx and ^ = VS{z) — Lz^- The continuous ECL is found as follows: 

= zjKzt 
= zJiH,-d,H,J 

= {Hjzt + Hlz,,) - {Hlz,, + zJd^H,^) 
= H,- (Hlz,),. 

The key step is the use of the product rule for the derivative of the flux Hj_^zt. Note 
that there is no requirement for K to be invertible. If K is invertible then one can 
solve and substitute for Zt to give the same flux as in Proposition 1. For the discrete 
gradient method we have 

= {z,- ZoVK{z, - zo)/At 
= {z,-zoV(H:-dJIZ) 

= {zi - ZqYH^ + {zi^ - zq^YH^^ - {{zi^ - zq-,Yh^^ + {zi - zoYd^H^J 
= H{zi, zi^) - H{zo, zox) - {H^J{zi - zo))x 
= AtH - (HZ^Atz),, 

establishing the result. □ 

We now consider the fully discrete case. The existence and derivation of con- 
servation laws for variational and Hamiltonian discretizations is the subject of dis- 
crete versions of Noether's theorem, see e.g. We restrict ourselves here to the 
observation that if the semidiscretization has a semidiscrete ECL, then a discrete 
gradient method applied to this semidiscretization will have a fully discrete ECL. In 
the following proposition, the extension to T-L = J H{z,Zx) is routine, we omit the 
^-dependence for clarity. 

Proposition 4. Let z{x,t) G M", let K be an n x n antisymmetric matrix and let 
y, be a Hamiltonian of the form H = J H{zx) dx with equations of motion 

bz 

Consider a finite difference discretization with = J H{zx) dx discretized to Hd = 
J2iH{Azi) where Azi := Zi^i — Zi. The semidiscretization 

Zt = KVHd 
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has a semidiscrete ECL and the discrete gradient methods applied to this semidis- 
cretization has a fully discrete ECL. 

Proof. The semidiscretization is 

= -KAVH{Azi-i) 

and the semidiscrete ECL is found as 

{H{Azi))t = VH{Az,y{Az,)t 

= -VHiAziY AKAVH{Azi^i) 

= -AVH{AziY KA\/H{Azi) - V H {Azif K AV H {Azi_i) 
= -A{VH{Azi)KAVH{Azi-i)) 

where the first term in the third fine is zero and the last hne follows because of the 
discrete product rule 

A(a,6,:_i) = 6jAaj + ajA6j_i. 
The parallel with the continuous ECL 

is clear. The discrete gradient method 

Zi — Zq 

At 

has fully discrete ECL 

H{zi) - H{zo) 



-KAH,^{Az,_,) 



= -KA{H,^{Azi) ' KAH,^{Az,_,)), 
the calculation proceeding exactly as in the spatially continuous case. □ 
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